Numerical investigation of the fractional diffusion wave equation with exponential kernel via cubic B-Spline approach

Splines are piecewise polynomials that are as smooth as they can be without forming a single polynomial. They are linked at specific points known as knots. Splines are useful for a variety of problems in numerical analysis and applied mathematics because they are simple to store and manipulate on a computer. These include, for example, numerical quadrature, function approximation, data fitting, etc. In this study, cubic B-spline (CBS) functions are used to numerically solve the time fractional diffusion wave equation (TFDWE) with Caputo-Fabrizio derivative. To discretize the spatial and temporal derivatives, CBS with θ-weighted scheme and the finite difference approach are utilized, respectively. Convergence analysis and stability of the presented method are analyzed. Some examples are used to validate the suggested scheme, and they show that it is feasible and fairly accurate.


Introduction
The area of mathematical analysis known as fractional calculus (FC) studies the applications of non-integer order integrals and derivatives.FC has generated significant interest in past decades.An extensive study on this topic is discussed in [1,2].Because of its wide uses in the numerous branches of engineering and science, the field of FC is highly regarded and significant.For instance, the FC has been implemented in the control theory [3], viscoelastic and viscoplastic flow [4], continuum mechanics [5], tumor development [6], transport problems [7], random walks [8], turbulence [9,10], coronavirus [11,12] and dynamical systems [13,14].Multitudinous authors [15][16][17][18][19] have called attention to the fractional integrals and derivatives that are certainly suitable for demonstrating memory and hereditary characteristics of several materials and methodologies that are handled by anomalous diffusion.In many cases, fractional problems characterised by fractional partial differential equations (FPDEs) behave more appropriately than their integer order equivalents.In general, it is impossible to obtain the exact result of the maximum of the FPDEs.As a result, many contemporary research articles focus on the quest for numerical approaches.In this paper, the following TFDWE with damping and reaction terms for different numerical outcomes shall be studied: @ g wðu; tÞ @t g þ d @wðu; tÞ @t þ swðu; tÞ À @ 2 wðu; tÞ @u 2 ¼ qðu; tÞ; g 2 ð1; 2Þ; u 2 ½a; b�; t 2 ½t 0 ; T�; ð1Þ with initial conditions (ICs): and boundary conditions (BCs): wða; tÞ ¼ � 1 ðtÞ; wðb; tÞ ¼ � 2 ðtÞ; ð3Þ where δ and σ are coefficients of the damping and reaction terms, respectively.q(u, t) is the source term and w t (u, t 0 ) represents the derivative of the function w(u, t) with respect to time t at t = t 0 .The Caputo-Fabrizio fractional derivative (CFFD) @ g @t g wðu; tÞ is described as: @ g @t g wðu; tÞ ¼ where R(γ) is the normalization function fulfils R(0) = R(1) = 1.Recently, several articles have discussed the CFFD's applications.For example, Cruz-Duarte et al. [20] propounded a closed formula for the Gaussian-based CFFD with signal processing applications.Shafiq et al. [21] demonstrated the numerical results of diffusion equation (DE) using spline functions.For solving optimal control problems, Mortezaee et al. [22] developed generalized fuzzy hyperbolic model.Zhou et al. [23] presented the model of solute transport and non-Darcian flow in porous media involving CFFD.The hepatitis E dynamics in CFFD sense has been presented by Khan et al. [24].Rahman et al. [25] studied the coronavirus (COVID-19) mathematical model with CFFD.Alshabanat et al. [26] presented the generalization of CFFD and its applications on electrical circuits.Recently, numerous numerical methodologies for resolving TFDWEs have been established.Ding and Li [27] developed two numerical techniques for solving the TFDWE with a reaction term.Avazzadeh et al. [28] used the radial basis functions to resolve the TFDWE.For fractional wave equation (FWE), Khader and Adel [29] introduced a Hermite formula based algorithm.Chatterjee et al. [30] proposed a method based on Bernstein polynomials for nonlinear TFDWE.Zeng [31] proposed two schemes for the TFDWE.The Legendre wavelets based computational method has been presented by Hooshmandasl et al. [32] to resolve the fractional sub-diffusion and TFDWEs.Ali et al. [33] applied implicit difference technique for computational solution of TFDWE.Sweilam et al. [34] investigated a computational study for variable order nonlinear FWE.Zhou and Xu [35] used Chebyshev wavelets collocation scheme for TFDWE.More classical work on TFDWE can be found in [36][37][38][39][40]. Koleva and Vulkov [41] used Jumarie's derivative to solve the Black-Scholes equation.Ganji et al. [42] proposed a new method for Klein-Gordon equation (KGE) with Caputo fractional derivative (CFD) through clique polynomials.Adomian decomposition technique has been utilized by Birajdar [43] for solving Navier-Stokes fractional equation.The numerical simulation for Fokker-Planck fractional equation has been presented by Mahdy [44].Hosseini et al. [45] proposed singular boundary technique for fractional DE with CFD.Akgu ¨l et al. [46] presented an accurate method for Lane-Emden fractional type equations.Babu et al. [47] developed Localized Differential Quadrature scheme for solving viscous Burgers' equation (BE).
Due to their simplicity, splines have been utilized by numerous researchers to resolve fractional differential equations.Wasim et al. [48] utilized hybrid B-spline scheme for solving Fisher and Huxley BEs.B-spline based computational methods have been proposed by Yaseen and Abbas [49,50] for fractional Burgers' and telegraph equations (TEs) with CFD.Akram et al. [51] presented extended CBS technique to solve non-linear TE involving CFD.Amin et al. [52] developed redefined extended CBS algorithm to resolve time-fractional KGE.Allen-Cahn equation by utilizing redefined CBS functions has been investigated by Khalid et al. [53] with CFD.Reaction-diffusion model with CFFD has been constructed by Akram et al. [54].Bsplines based technique for time-fractional DE with CFD has been presented by Yaseen et al. [55].Akram et al. [56] developed a numerical approach which depends on extended CBS function for Fisher equation with CFD.Shafiq et al. [57] resolved the Burgers' equations involving ABFD numerical using CBS.Kamran et al. [58] presented numerical simulation for BBM-BE with CFD using CBS functions.Shafiq et al. presented the CBS based algorithm to resolve the TFADE in the article [59].
The proposed attempt has been inspired by recent developments in the investigation of computational solutions to the TFDWE.As a result of its extensive use, our study's objective is to apply CBS to the TFDWE.In the past, many researchers utilized B-spline techniques to resolve TFDWE, but no one has ever employed CBS approximations with CFFD.To solve the TFDWE, the CBS with θ-weighted technique is used.Analyses of the scheme's convergence and stability are also carried out.By giving answers to a few numerical issues, the propounded approach's efficacy and applicability are illustrated.We find that our proposed technique yields effective results by contrasting the attained numerical solutions with the analytical solutions.The author is aware of no studies that have been conducted on the suggested algorithm for TFDWE.
The outline of this research paper are: Parseval's identity (PI) and CBS functions are included in section 2, and the proposed algorithm is discussed in section 3. The stability of current methodology and analysis of convergence are covered in sections 4 and 5, respectively.In section 6, the effectiveness and validity of the suggested technique are investigated, and finally, section 7 summarizes the conclusion.

Preliminaries
Definition 1.If g 2 L 2 ½a; b�, then PI is defined as [60]: where ĝ ðmÞ ¼ R b a g ðuÞe 2pimu du is Fourier transform for all integers m.

Basis functions
Consider a = u 0 < u 1 < � � � < u N = b be the equivalent partition of [a, b] at the knots u r = u 0 + rh, r = 0, 1, � � �, N, where h ¼ bÀ a N .The CBS functions are expressed as [57]: where u is the variable.The CBS has geometrical characteristics like non-negativity, convex hull property, symmetry, partition of unity and geometric invariability [53].Additionally, S −1 , S 0 , � � �, S N+1 have been constructed.For w(u, t), the approximation W(u, t) in terms of CBS can be assumed as [59]: where control points ϑ r (t) to be calculated at every temporal stage.Eqs ( 6) and ( 7) provide the following approximations at the nodal points: 3 Description of numerical scheme Using forward difference approach, Eq (9) becomes where m ¼ 1 À exp ðÀ g 2À g DtÞ and k s ¼ exp ðÀ g 2À g sDtÞ.Moreover, the truncation error F nþ1 Dt is given as where Ϝ is a constant.It is straightforward to check that Using (10) and θ-weighted scheme, Eq (1) becomes mRðgÞ gðDtÞ will appear when n = s or n = 0. To remove w −1 , we use the IC to attain Employing the CBS approximation and its necessary derivatives at knot u r in (13), then Substituting ( 8) in ( 15), we get This system (16) has N + 1 linear equations involving N + 3 unknowns.To attain a consistent system, two additional equations are acquired by utilizing the BCs (3).Hence, a matrix system of dimension (N + 3) × (N + 3) is attained which can be uniquely solved utilizing any appropriate algorithm.Before using (16), the initial vector achieved through the ICs as The matrix form of Eq ( 17) is expressed as where and

:
Any numerical algorithm can be used to solve Eq (18) for ϑ 0 .Mathematica is utilized to accomplish the numerical results.

The stability of presented scheme
If the error does not increase while the computation is in progress, it can be inferred that the numerical method is stable [61].Fourier method is employed for the stability of presented scheme.For this, suppose that B n r and Bn r represent the growth factors analytically and numerically, respectively.The error φ n r can be written as Thus from Eq (16), we obtain From ICs and BCs, we can write and The grid function is stated as: The φ n (u) in the Fourier mode can be presented as: where Implementing k.k 2 norm, we attain hjφ n r j Using Parseval's identity (5), we achieve Hence, we acquire Assume that the Eqs ( 19)-( 21) possess the solution in Fourier sense as: and ρ is any real number.Substituting (26) in (19) and simplifying, we achieve Utilizing the relation e iρh + e −iρh = 2cos(ρh) and simplifying, we get where @ ¼ % %þd 0 and ε ¼ 1 þ 2 cos ðrhÞðsh 2 À 6Þþ4sh 2 þ12 ð%h 2 þd 0 h 2 Þð2 cos ðrhÞþ4Þ .Obviously ε � 1. Proposition 4.1.If χ n is the solution for Eq (28), Proof.To prove this, mathematical induction is utilized.For n = 0, Eq (28) gives The proposed method ( 16) is unconditionally stable.Proof.Employing Eq (25) and Proposition 4.1, we get Thus, the proposed computational technique is stable unconditionally.
Proof.Let W ðu; tÞ ¼ P Nþ1 r¼À 1 x n r ðtÞS r ðuÞ be the estimated CBS for W(u, t).Through the triangular inequality, we get kwðu; tÞ À Wðu; tÞk 1 � kwðu; tÞ À W ðu; tÞk 1 þ k W ðu; tÞ À Wðu; tÞk 1 : With the aid of Theorem 2 for r = 0, we attain The present method possesses collocation conditions Lw(u r , t) = LW(u r , t) = q(u r , t), r = 0, 1, � � �, N. Consider L W ðu r ; tÞ ¼ qðu r ; tÞ: Therefore, the difference equation Lð W ðu r ; tÞ À Wðu r ; tÞÞ can be stated at time level n as The BCs can be given as where From inequality (29), we achieve jz n r j, e n r ¼ jO n r j and e n ¼ max 0�r�N je n r j.When n = 0 and using the expression (14), Eq (33) becomes where r = 0, 1, � � �, N. Using IC, e 0 = 0: Taking absolute values of O 1 r , z 1 r and appropriately small h, we have We get the values of e 1 À 1 and e 1 Nþ1 through BCs: which implies where ξ 1 is not depending on h.Now, mathematical induction is utilized for the proof of this theorem.Suppose that e y r � x y h 2 is true for 1 � y � n and ξ = max{ξ y : y = 0, 1, � � �, n}, then from Eq (33), we acquire Again, taking absolute values of O nþ1 r and z nþ1 r , we obtain Similarly, we get the values of e nþ1 À 1 and e nþ1 Nþ1 from the boundary conditions Hence, for all n, we acquire x. Theorem 4. The TFDWE is convergent with ICs and BCs.Proof.Assume that w(u, t) and W(u, t) are analytical and approximate results for TFDWE, respectively.Therefore, the preceding theorem and relation (11) justify that there exist constants x and Ϝ such that Thus, the proposed method is second order convergent.

Numerical results and discussion
In this section, numerical results of some experiments using proposed approach are demonstrated.To examine the validity of the proposed scheme, we employ error norms L 2 and L 1 as  where qðu; tÞ ¼ 2 RðgÞ The w(u, t) = (t 2 − t)sin(πu) is the exact solution.The approximate results and absolute error for distinct u values at t = 0.2 for Example 6.1 are demonstrated in Table 1.For various time stages, error norms L 2 and L 1 are displayed in Table 2. Table 3 contains the analysis of error norms and convergence order.Tables 4 and 5 Example 6.2.Consider the TFDWE @ g wðu; tÞ @t g þ @wðu; tÞ @t À @ 2 wðu; tÞ @u The w(u, t) = t 2 u(1−u) is the exact solution.
where qðu; tÞ The w(u, t) = t 2 sinh(u) is the exact solution.Table 10 displays computational outcomes and absolute error of Example 6.3 at various spatial grid values.Table 11 expounds the analysis of error norms and convergence order in the spatial direction.Table 12 describes the error norm at several values of Δt.Table 13 depicts the error norm for γ = 1.5, 1.9 at numerous time stages.

Conclusion
In this paper, we addressed the problem of finding numerical solutions to time FPDE.For this purpose, B-splines were utilized to build up a collocation technique for TFDWE.The procedure employed the typical FDM to approximate the time fractional derivative, whereas the derivative in space was discretized utilizing the cubic B-splines.Further, the effective numerical scheme to examine the numerical solutions of TFDWE involving damping and reaction terms has been presented.We have used the CBS functions and θ-weighted scheme with CFFD.The proposed scheme is stable unconditionally possessing second order temporal and spatial convergence.Three numerical problems have been analyzed.Numerical and graphical comparison uncovers that the provided method is accurate and computationally very effective.
In future, we may consider the higher dimensional and variable order FPDEs for numerical solutions via spline functions.

Fig 3 .
Fig 3. 2D and 3D error profiles for Δt = 0.01, N = 180,γ = 1.5, u 2 [0, 1] and t = 1 for Example 6.1.(a) 2D error function (b) 3D error function.https://doi.org/10.1371/journal.pone.0295525.g003 Fig 4 exhibits the relation between analytical solutions and numerical outcomes at different temporal directions.3D precision of the existing approach is displayed by graphs of numerical results and analytic solutions in Fig 5. Fig 6 demonstrates the 2D and 3D errors description.Tables and figures dipict that the presented numerical results are adaptable with exact solutions.The >

Table 4 . Error norm at t = 0.4 and γ = 1.7 for Example 6.1.
Table 6 displays the absolute error and the approximate results of Example 6.2 at various spatial grid values.At numerous time stages, error norms L 2 and L 1 are represented in Table 7. Table 8 comprises the analysis of error norms and convergence order.Table 9 describes the error norm at several values of Δt.